# Load Necessary Packages
library(foreign)
library(stargazer)
library(Amelia)
library(Zelig)
library(dplyr)
library(MuMIn)
library(haven)
library(ggplot2)
library(grid)
library(gridExtra)
library(psych)
library(memisc)
library(ggpubr)

# Rescale Command
rescale <- function(x) {
  if(min(x,na.rm=T)>=0){
    x <- (x-min(x,na.rm=T))/(max(x,na.rm=T) - min(x,na.rm=T))
  } else{
    x <- x + abs(0 - min(x,na.rm=T))
    x <- (x-min(x,na.rm=T))/(max(x,na.rm=T) - min(x,na.rm=T))
  }
  print(x)
}

# Load Data
setwd("/users/josephphillips/Dropbox/Projects/Sorting and Affective Polarization/")
anes.cumulative <- read.dta("anes_timeseries_cdf.dta",convert.factors=T)
anes.cumulative <- subset(anes.cumulative,VCF0303=="1. Democrats (including leaners)" | VCF0303=="3. Republicans (including leaners)")

# Clean Feeling Thermometers
anes.cumulative$VCF0201[anes.cumulative$VCF0201>97] <- NA
anes.cumulative$VCF0202[anes.cumulative$VCF0202>97] <- NA
anes.cumulative$VCF0101[anes.cumulative$VCF0101==0] <- NA
anes.cumulative$VCF0316[anes.cumulative$VCF0316>5] <- NA
anes.cumulative$VCF0320[anes.cumulative$VCF0320>5] <- NA

# Create Affective Polarization Variable
anes.cumulative$polariz.old <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0201-anes.cumulative$VCF0202),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0202-anes.cumulative$VCF0201),NA))
anes.cumulative$polariz.new <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0218-anes.cumulative$VCF0224),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0224-anes.cumulative$VCF0218),NA))
anes.cumulative$polariz.likes <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0316-anes.cumulative$VCF0320),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0320-anes.cumulative$VCF0316),NA))
anes.cumulative$polariz.combined <- ifelse(anes.cumulative$VCF0004<1978,anes.cumulative$polariz.old,anes.cumulative$polariz.new)
anes.cumulative$inparty.new <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0218),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0224),NA))
anes.cumulative$outparty.new <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0224),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0218),NA))
anes.cumulative$inparty.old <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0201),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0202),NA))
anes.cumulative$outparty.old <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0202),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0201),NA))
anes.cumulative$inparty.likes <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0316),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0320),NA))
anes.cumulative$outparty.likes <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0320),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0316),NA))
anes.cumulative$inparty.combined <- ifelse(anes.cumulative$VCF0004<1978,anes.cumulative$inparty.old,anes.cumulative$inparty.new)
anes.cumulative$outparty.combined <- ifelse(anes.cumulative$VCF0004<1978,anes.cumulative$outparty.old,anes.cumulative$outparty.new)

# Predictors of Affective Polarization
## PID
anes.cumulative$pid <- ifelse(anes.cumulative$VCF0301=="0. DK; NA; other; refused to answer; no Pre IW",NA,as.numeric(anes.cumulative$VCF0301) - 1)
## Strong Identifier
anes.cumulative$strong <- ifelse(anes.cumulative$pid==1 | anes.cumulative$pid==7,1,0)
## Age
anes.cumulative$age <- anes.cumulative$VCF0101-18
anes.cumulative$agesq <- anes.cumulative$age * anes.cumulative$age
## Gender
anes.cumulative$gender <- ifelse(anes.cumulative$VCF0104==0,NA,ifelse(anes.cumulative$VCF0104==8,0.5,ifelse(anes.cumulative$VCF0104==2,1,0)))
## Race
anes.cumulative$black <- ifelse(anes.cumulative$VCF0106=="0. Missing, pre-1966 data" | anes.cumulative$VCF0106=="9. Missing, DK/REF/NA",NA,ifelse(anes.cumulative$VCF0106=="2. Black non-Hispanic",1,0))
anes.cumulative$blacksort <- ifelse(anes.cumulative$pid<4 & anes.cumulative$black==1,1,ifelse(anes.cumulative$pid>4 & anes.cumulative$black==1,-1,0))
## Education
anes.cumulative$education <- ifelse(anes.cumulative$VCF0140=="8. DK" | anes.cumulative$VCF0140=="9. NA; RF; no Pre IW; short-form 'new' Cross Section",NA,as.numeric(anes.cumulative$VCF0140))
## South
anes.cumulative$south <- ifelse(anes.cumulative$VCF0112=="3. South (AL, AR, DE, D.C., FL, GA, KY, LA, MD, MS, NC",1,0)
## Ideological Sorting
anes.cumulative$ideo <- ifelse(anes.cumulative$VCF0803=="0. NA; no Post IW; form III,IV (1972); R not" | anes.cumulative$VCF0803=="9. DK; haven't thought much about it",NA,as.numeric(anes.cumulative$VCF0803)-1)
anes.cumulative$stronglib <- ifelse(anes.cumulative$VCF0803=="9. DK; haven't thought much about it" | anes.cumulative$VCF0803=="0. NA; no Post IW; form III,IV (1972); R not",NA,ifelse(anes.cumulative$VCF0803=="1. Extremely liberal",1,0))
anes.cumulative$strongcon <- ifelse(anes.cumulative$VCF0803=="9. DK; haven't thought much about it" | anes.cumulative$VCF0803=="0. NA; no Post IW; form III,IV (1972); R not",NA,ifelse(anes.cumulative$VCF0803=="7. Extremely conservative",1,0))
anes.cumulative$strongsort <- ifelse(anes.cumulative$pid<4 & anes.cumulative$stronglib==1,1,ifelse(anes.cumulative$pid<4 & anes.cumulative$strongcon==1,-1,ifelse(anes.cumulative$pid>4 & anes.cumulative$strongcon==1,1,ifelse(anes.cumulative$pid>4 & anes.cumulative$stronglib==1,-1,0))))
anes.cumulative$ideosort <- rescale(-abs(anes.cumulative$pid-anes.cumulative$ideo))
## Church Attendance
anes.cumulative$attendance <- as.numeric(anes.cumulative$VCF0130)
anes.cumulative$attendance[anes.cumulative$attendance==1 | anes.cumulative$attendance==6 | anes.cumulative$attendance==7] <- 1
anes.cumulative$attendance[anes.cumulative$attendance==5] <- -2
anes.cumulative$attendance[anes.cumulative$attendance==4] <- -3
anes.cumulative$attendance[anes.cumulative$attendance==3] <- -4
anes.cumulative$attendance[anes.cumulative$attendance==2] <- -5
anes.cumulative$attendance[anes.cumulative$attendance==-2] <- 2
anes.cumulative$attendance[anes.cumulative$attendance==-3] <- 3
anes.cumulative$attendance[anes.cumulative$attendance==-4] <- 4
anes.cumulative$attendance[anes.cumulative$attendance==-5] <- 5
anes.cumulative$attendance[anes.cumulative$attendance>7] <- NA
anes.cumulative$attendance <- (rescale(anes.cumulative$attendance)*2)-1
anes.cumulative$churchsort <- ifelse(anes.cumulative$pid>4,anes.cumulative$attendance,-anes.cumulative$attendance)
## Social Sorting
anes.cumulative$socialsort <- rescale(anes.cumulative$strong+anes.cumulative$blacksort+anes.cumulative$strongsort+anes.cumulative$churchsort)
## Interest (Media Proxy)
anes.cumulative$media <- ifelse(anes.cumulative$VCF0310=="0. NA; no Pre IW; version B (2008)" | anes.cumulative$VCF0310=="9. DK",NA,as.numeric(anes.cumulative$VCF0310)-2)
## Year
anes.cumulative$year <- anes.cumulative$VCF0004
anes.cumulative$cable <- ifelse(anes.cumulative$year>=1996,1,0)
anes.cumulative$year <- as.factor(anes.cumulative$VCF0004)

## APC Variables
anes.cumulative$agegroup <- ifelse(anes.cumulative$age<=6,20,ifelse(anes.cumulative$age>6 & anes.cumulative$age<=11,25,ifelse(anes.cumulative$age>11 & anes.cumulative$age<=16,30,ifelse(anes.cumulative$age>16 & anes.cumulative$age<=21,35,ifelse(anes.cumulative$age>21 & anes.cumulative$age<=26,40,ifelse(anes.cumulative$age>26 & anes.cumulative$age<=31,45,ifelse(anes.cumulative$age>31 & anes.cumulative$age<=36,50,ifelse(anes.cumulative$age>36 & anes.cumulative$age<=41,55,ifelse(anes.cumulative$age>41 & anes.cumulative$age<=46,60,ifelse(anes.cumulative$age>46 & anes.cumulative$age<=51,65,ifelse(anes.cumulative$age>51 & anes.cumulative$age<=56,70,ifelse(anes.cumulative$age>56 & anes.cumulative$age<=61,75,ifelse(anes.cumulative$age>61 & anes.cumulative$age<=66,80,ifelse(anes.cumulative$age>66 & anes.cumulative$age<=71,85,ifelse(anes.cumulative$age>71,90,NA)))))))))))))))
anes.cumulative$yeargroup <- ifelse(anes.cumulative$VCF0004>1975 & anes.cumulative$VCF0004<=1980,1980,ifelse(anes.cumulative$VCF0004>1980 & anes.cumulative$VCF0004<=1985,1985,ifelse(anes.cumulative$VCF0004>1985 & anes.cumulative$VCF0004<=1990,1990,ifelse(anes.cumulative$VCF0004>1990 & anes.cumulative$VCF0004<=1995,1995,ifelse(anes.cumulative$VCF0004>1995 & anes.cumulative$VCF0004<=2000,2000,ifelse(anes.cumulative$VCF0004>2000 & anes.cumulative$VCF0004<=2005,2005,ifelse(anes.cumulative$VCF0004>2005 & anes.cumulative$VCF0004<=2010,2010,ifelse(anes.cumulative$VCF0004>2010 & anes.cumulative$VCF0004<=2015,2015,ifelse(anes.cumulative$VCF0004>2015,2020,NA)))))))))
anes.cumulative$yeargroup.old <- ifelse(anes.cumulative$VCF0004>1960 & anes.cumulative$VCF0004<=1965,1965,ifelse(anes.cumulative$VCF0004>1965 & anes.cumulative$VCF0004<=1970,1970,ifelse(anes.cumulative$VCF0004>1970 & anes.cumulative$VCF0004<=1975,1975,ifelse(anes.cumulative$VCF0004>1975 & anes.cumulative$VCF0004<=1980,1980,ifelse(anes.cumulative$VCF0004>1980 & anes.cumulative$VCF0004<=1985,1985,ifelse(anes.cumulative$VCF0004>1985 & anes.cumulative$VCF0004<=1990,1990,ifelse(anes.cumulative$VCF0004>1990 & anes.cumulative$VCF0004<=1995,1995,ifelse(anes.cumulative$VCF0004>1995 & anes.cumulative$VCF0004<=2000,2000,ifelse(anes.cumulative$VCF0004>2000 & anes.cumulative$VCF0004<=2005,2005,ifelse(anes.cumulative$VCF0004>2005 & anes.cumulative$VCF0004<=2010,2010,ifelse(anes.cumulative$VCF0004>2010 & anes.cumulative$VCF0004<=2015,2015,ifelse(anes.cumulative$VCF0004>2015,2020,NA))))))))))))
anes.cumulative$yeargroup.likes <- ifelse(anes.cumulative$VCF0004>1950 & anes.cumulative$VCF0004<=1955,1955,ifelse(anes.cumulative$VCF0004>1955 & anes.cumulative$VCF0004<=1960,1960,ifelse(anes.cumulative$VCF0004>1960 & anes.cumulative$VCF0004<=1965,1965,ifelse(anes.cumulative$VCF0004>1965 & anes.cumulative$VCF0004<=1970,1970,ifelse(anes.cumulative$VCF0004>1970 & anes.cumulative$VCF0004<=1975,1975,ifelse(anes.cumulative$VCF0004>1975 & anes.cumulative$VCF0004<=1980,1980,ifelse(anes.cumulative$VCF0004>1980 & anes.cumulative$VCF0004<=1985,1985,ifelse(anes.cumulative$VCF0004>1985 & anes.cumulative$VCF0004<=1990,1990,ifelse(anes.cumulative$VCF0004>1990 & anes.cumulative$VCF0004<=1995,1995,ifelse(anes.cumulative$VCF0004>1995 & anes.cumulative$VCF0004<=2000,2000,ifelse(anes.cumulative$VCF0004>2000 & anes.cumulative$VCF0004<=2005,2005,ifelse(anes.cumulative$VCF0004>2005 & anes.cumulative$VCF0004<=2010,2010,ifelse(anes.cumulative$VCF0004>2010 & anes.cumulative$VCF0004<=2015,2015,ifelse(anes.cumulative$VCF0004>2015,2020,NA))))))))))))))
anes.cumulative$yeargroup.old.na <- ifelse(anes.cumulative$yeargroup.old>1985,NA,anes.cumulative$yeargroup.old)
anes.cumulative$yeargroup.likes.na <- ifelse(anes.cumulative$yeargroup.likes>2005,NA,anes.cumulative$yeargroup.likes)
anes.cumulative$cohort <- anes.cumulative$yeargroup-anes.cumulative$agegroup
anes.cumulative$cohort.old <- anes.cumulative$yeargroup.old.na-anes.cumulative$agegroup
anes.cumulative$cohort.combined <- anes.cumulative$yeargroup.old-anes.cumulative$agegroup
anes.cumulative$cohort.likes <- anes.cumulative$yeargroup.likes.na-anes.cumulative$agegroup
anes.cumulative$cohort.likes.full <- anes.cumulative$yeargroup.likes-anes.cumulative$agegroup
anes.cumulative$agegroup <- as.factor(anes.cumulative$agegroup)
anes.cumulative$yeargroup <- as.factor(anes.cumulative$yeargroup)
anes.cumulative$yeargroup.old <- as.factor(anes.cumulative$yeargroup.old)
anes.cumulative$yeargroup.likes <- as.factor(anes.cumulative$yeargroup.likes)
anes.cumulative$yeargroup.old.na <- as.factor(anes.cumulative$yeargroup.old.na)
anes.cumulative$yeargroup.likes.na <- as.factor(anes.cumulative$yeargroup.likes.na)
anes.cumulative$cohort <- as.factor(anes.cumulative$cohort)
anes.cumulative$cohort.old <- as.factor(anes.cumulative$cohort.old)
anes.cumulative$cohort.combined <- as.factor(anes.cumulative$cohort.combined)
anes.cumulative$cohort.likes <- as.factor(anes.cumulative$cohort.likes)

anes.cumulative.full <- anes.cumulative
anes.cumulative.dem <- subset(anes.cumulative.full,pid<=3)
anes.cumulative.rep <- subset(anes.cumulative.full,pid>=5)

# Affective Polarization
## Raw Model, All Partisans
new.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
summary(new.all.raw$model) ## Age and period effects
new.all.raw.cohort <- APCI::cohortdeviation(new.all.raw$A,new.all.raw$P,new.all.raw$C,model=new.all.raw$model,covariate=c())
new.all.raw.cohort$cohort_average ## cohort means
new.all.raw.cohort$cohort_slope ## cohort slopes
## Model Controlling for Partisan Strength
newpid.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong'),family='gaussian')
summary(newpid.all.raw$model) # Age effects
## Raw Model, Democrats
new.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
## Raw Model, Republicans
repnot90 <- subset(anes.cumulative.rep,agegroup!="90")
repnot90$agegroup <- as.factor(ifelse(repnot90$agegroup=="20","20",ifelse(repnot90$agegroup=="25","25",ifelse(repnot90$agegroup=="30","30",ifelse(repnot90$agegroup=="35","35",ifelse(repnot90$agegroup=="40","40",ifelse(repnot90$agegroup=="45","45",ifelse(repnot90$agegroup=="50","50",ifelse(repnot90$agegroup=="55","55",ifelse(repnot90$agegroup=="60","60",ifelse(repnot90$agegroup=="65","65",ifelse(repnot90$agegroup=="70","70",ifelse(repnot90$agegroup=="75","75",ifelse(repnot90$agegroup=="80","80","85"))))))))))))))
new.rep.raw <- APCI::temp_model(data=repnot90,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
## Model Controlling for Sorting, Democrats
newsort.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
## Model Controlling for Sorting, Republicans
newsort.rep.raw <- APCI::temp_model(data=repnot90,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
## Figure 1 (Age Effects, IE coefficients and standard errors lifted from .do file)
figure1data <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                          effect_apci=c(NA,summary(new.all.raw$model)$coefficients[2:15,1]),
                          se_apci=c(NA,summary(new.all.raw$model)$coefficients[2:15,2]),
                          effect_ie=c(-6.415,-5.220,-5.128,-3.457,-1.722,0.773,-1.505,1.156,1.686,3.211,4.113,2.281,2.644,2.618,6.510),
                          se_ie=c(0.953,0.844,0.751,0.682,0.663,0.653,0.640,0.658,0.715,0.779,0.906,1.057,1.268,1.723,3.083))
figure1data$lci_apci <- figure1data$effect_apci - (1.96*figure1data$se_apci)
figure1data$uci_apci <- figure1data$effect_apci + (1.96*figure1data$se_apci)
figure1data$lci_ie <- figure1data$effect_ie - (1.96*figure1data$se_ie)
figure1data$uci_ie <- figure1data$effect_ie + (1.96*figure1data$se_ie)
figure1data$color_ie <- ifelse(figure1data$lci_ie<0 & figure1data$uci_ie>0,"Gray50","Black")
figure1data$color_apci <- ifelse(figure1data$lci_apci<0 & figure1data$uci_apci>0,"Gray50","Black")
figure1data$group <- 1

figure1.a <- ggplot(figure1data,aes(x=age,y=effect_apci,color=color_apci,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_apci,ymax=uci_apci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + ylab("Deviation from Grand Mean") + xlab("Age")
figure1.b <- ggplot(figure1data,aes(x=age,y=effect_ie,color=color_ie,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_ie,ymax=uci_ie),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + ylab("Deviation from Grand Mean") + xlab("Age")
figure1 <- ggarrange(figure1.a,figure1.b,nrow=1)
## Figure 2 (Period Effects, IE coefficients and standard errors lifted from .do file)
figure2data <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                          effect_ie=c(NA,summary(new.all.raw$model)$coefficients[16:23,1]),
                          se_ie=c(NA,summary(new.all.raw$model)$coefficients[16:23,2]),
                          effect_apci=c(NA,-7.059,-2.858,-4.962,-4.132,-2.413,2.393,2.836,8.915),
                          se_apci=c(NA,1.045,0.687,0.795,1.145,0.741,1.006,1.021,1.350))
figure2data$lci_ie <- figure2data$effect_ie - (1.96*figure2data$se_ie)
figure2data$uci_ie <- figure2data$effect_ie + (1.96*figure2data$se_ie)
figure2data$lci_apci <- figure2data$effect_apci - (1.96*figure2data$se_apci)
figure2data$uci_apci <- figure2data$effect_apci + (1.96*figure2data$se_apci)
figure2data$color_ie <- ifelse(figure2data$lci_ie<0 & figure2data$uci_ie>0,"Gray50","Black")
figure2data$color_apci <- ifelse(figure2data$lci_apci<0 & figure2data$uci_apci>0,"Gray50","Black")
figure2data$group <- 1

figure2.a <- ggplot(figure2data,aes(x=period,y=effect_apci,color=color_apci,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_apci,ymax=uci_apci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + ylab("Deviation from Grand Mean") + xlab("Year")
figure2.b <- ggplot(figure2data,aes(x=period,y=effect_ie,color=color_ie,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_ie,ymax=uci_ie),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + ylab("Deviation from Grand Mean") + xlab("Year")
ggarrange(figure2.a,figure2.b,nrow=1)
## Figure 4, Left Panel (Cohort Means, IE analyses lifted from .do file)
figure3data <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                          effect_ie=c(20.142,-5.557,-2.674,0.827,-0.447,-3.783,-1.634,-3.493,-2.945,-1.647,-3.417,-2.301,-1.527,-0.900,0.542,0.005,-0.648,-0.919,1.143,1.340,0.558,2.752,-0.767),
                          se_ie=c(9.335,5.730,3.190,2.530,1.947,1.689,1.531,1.384,1.240,1.161,1.041,0.904,0.772,0.664,0.577,0.563,0.625,0.746,0.860,0.991,1.077,1.238,2.001),
                          effect_apci=c(as.numeric(new.all.raw.cohort$cohort_average$cohort_average)),
                          se_apci=c(as.numeric(new.all.raw.cohort$cohort_average$cohort_average_se)))
figure3data$lci_ie <- figure3data$effect_ie - (1.96*figure3data$se_ie)
figure3data$uci_ie <- figure3data$effect_ie + (1.96*figure3data$se_ie)
figure3data$lci_apci <- figure3data$effect_apci - (1.96*figure3data$se_apci)
figure3data$uci_apci <- figure3data$effect_apci + (1.96*figure3data$se_apci)
figure3data$color_ie <- ifelse(figure3data$lci_ie<0 & figure3data$uci_ie>0,"Gray50","Black")
figure3data$color_apci <- ifelse(figure3data$lci_apci<0 & figure3data$uci_apci>0,"Gray50","Black")
figure3data$group <- 1

figure3.a <- ggplot(figure3data,aes(x=cohort,y=effect_apci,color=color_apci,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_apci,ymax=uci_apci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + ylab("Deviation from Grand Mean") + xlab("Birth Cohort")
figure3.b <- ggplot(figure3data,aes(x=cohort,y=effect_ie,color=color_ie,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_ie,ymax=uci_ie),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + ylab("Deviation from Grand Mean") + xlab("Birth Cohort")
ggarrange(figure3.a,figure3.b,nrow=1)
## Figure 5 (Cohort Slopes)
figure5data <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                          effect=c(as.numeric(new.all.raw.cohort$cohort_slope$cohort_slope)),
                          se=c(as.numeric(new.all.raw.cohort$cohort_slope$cohort_slope_se)))
figure5data$lci <- figure5data$effect - (1.96*figure5data$se)
figure5data$uci <- figure5data$effect + (1.96*figure5data$se)
figure5data$color <- ifelse(figure5data$lci<0 & figure5data$uci>0,"Gray50","Black")
figure5data$group <- 1

ggplot(figure5data,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + xlab("Birth Cohort") + ylab("Deviation from Grand Mean")
## Figure 9 (Age Effects without and with controlling for PID Strength)
figure9data <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                          effect_raw=c(summary(new.all.raw$model)$coefficients[2:15,1]),
                          se_raw=c(summary(new.all.raw$model)$coefficients[2:15,2]),
                          effect_controlled=c(summary(newpid.all.raw$model)$coefficients[3:16,1]),
                          se_controlled=c(summary(newpid.all.raw$model)$coefficients[3:16,2]))
figure9data$lci_raw <- figure9data$effect_raw - (1.96*figure9data$se_raw)
figure9data$uci_raw <- figure9data$effect_raw + (1.96*figure9data$se_raw)
figure9data$lci_controlled <- figure9data$effect_controlled - (1.96*figure9data$se_controlled)
figure9data$uci_controlled <- figure9data$effect_controlled + (1.96*figure9data$se_controlled)
figure9data$color_raw <- ifelse(figure9data$lci_raw<0 & figure9data$uci_raw>0,"Gray50","Black")
figure9data$color_controlled <- ifelse(figure9data$lci_controlled<0 & figure9data$uci_controlled>0,"Gray50","Black")
figure9data$group <- 1

figure9.a <- ggplot(figure9data,aes(x=age,y=effect_raw,color=color_raw,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_raw,ymax=uci_raw),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Raw") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + xlab("Age") + ylab("Deviation from Grand Mean")
figure9.b <- ggplot(figure9data,aes(x=age,y=effect_controlled,color=color_controlled,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_controlled,ymax=uci_controlled),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Controlling") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + xlab("Age") + ylab("Deviation from Grand Mean")
ggarrange(figure9.a,figure9.b,nrow=1)
## Figure 10 (Period Effects with and without controlling for Sorting, Democrats)
figure10data <- data.frame(period=c(1980,1985,1990,1995,2000,2005,2010,2015),
                           effect_raw=c(summary(new.dem.raw$model)$coefficients[16:23,1]),
                           se_raw=c(summary(new.dem.raw$model)$coefficients[16:23,2]),
                           effect_controlled=c(summary(newsort.dem.raw$model)$coefficients[18:25,1]),
                           se_controlled=c(summary(newsort.dem.raw$model)$coefficients[18:25,2]))

figure10data$lci_raw <- figure10data$effect_raw - (1.96*figure10data$se_raw)
figure10data$uci_raw <- figure10data$effect_raw + (1.96*figure10data$se_raw)
figure10data$lci_controlled <- figure10data$effect_controlled - (1.96*figure10data$se_controlled)
figure10data$uci_controlled <- figure10data$effect_controlled + (1.96*figure10data$se_controlled)
figure10data$color_raw <- ifelse(figure10data$lci_raw<0 & figure10data$uci_raw>0,"Gray50","Black")
figure10data$color_controlled <- ifelse(figure10data$lci_controlled<0 & figure10data$uci_controlled>0,"Gray50","Black")
figure10data$group <- 1

figure10.a <- ggplot(figure10data,aes(x=period,y=effect_raw,color=color_raw,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_raw,ymax=uci_raw),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Raw") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + xlab("Year") + ylab("Deviation from Grand Mean")
figure10.b <- ggplot(figure10data,aes(x=period,y=effect_controlled,color=color_controlled,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_controlled,ymax=uci_controlled),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Controlling") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + xlab("Year") + ylab("Deviation from Grand Mean")
figure10 <- ggarrange(figure10.a,figure10.b,nrow=1)
## Figure 11 (Period Effects with and without controlling for sorting, Republicans)
figure11data <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85),
                           effect_raw=c(summary(new.rep.raw$model)$coefficients[2:14,1]),
                           se_raw=c(summary(new.rep.raw$model)$coefficients[2:14,2]),
                           effect_controlled=c(summary(newsort.rep.raw$model)$coefficients[4:16,1]),
                           se_controlled=c(summary(newsort.rep.raw$model)$coefficients[4:16,2]))

figure11data$lci_raw <- figure11data$effect_raw - (1.96*figure11data$se_raw)
figure11data$uci_raw <- figure11data$effect_raw + (1.96*figure11data$se_raw)
figure11data$lci_controlled <- figure11data$effect_controlled - (1.96*figure11data$se_controlled)
figure11data$uci_controlled <- figure11data$effect_controlled + (1.96*figure11data$se_controlled)
figure11data$color_raw <- ifelse(figure11data$lci_raw<0 & figure11data$uci_raw>0,"Gray50","Black")
figure11data$color_controlled <- ifelse(figure11data$lci_controlled<0 & figure11data$uci_controlled>0,"Gray50","Black")
figure11data$group <- 1

figure11.a <- ggplot(figure11data,aes(x=age,y=effect_raw,color=color_raw,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_raw,ymax=uci_raw),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Raw") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + xlab("Year") + ylab("Deviation from Grand Mean")
figure11.b <- ggplot(figure11data,aes(x=age,y=effect_controlled,color=color_controlled,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci_controlled,ymax=uci_controlled),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Controlling") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylim(-22,44) + xlab("Year") + ylab("Deviation from Grand Mean")
figure11 <- ggarrange(figure11.a,figure11.b,nrow=1)

# Comparing Affective Polarization in ANES and YPSP
## Import YPSP and code affective polarization
setwd("/Users/josephphillips/Dropbox/Projects/Dissertation/REAL Chapter 1 Stability/")
ypsp <- as.data.set(spss.portable.file('Jennings Niemi.por'))
ypsp$pid73 <- ifelse(as.numeric(ypsp$V482)<7,as.numeric(ypsp$V482)+1,NA) # PID in 1973
ypsp$pid82 <- ifelse(as.numeric(ypsp$V1608)<8,as.numeric(ypsp$V1608)+1,NA) # PID in 1982
ypsp$pid97 <- ifelse(as.numeric(ypsp$V5754)<7,as.numeric(ypsp$V5754)+1,NA) # PID in 1997
ypsp$dem73 <- as.numeric(ypsp$V432) # Democrats FT, 1973
ypsp$rep73 <- as.numeric(ypsp$V440) # Republicans FT, 1973
ypsp$dem82 <- as.numeric(ypsp$V1520) # Democrats FT, 1982
ypsp$rep82 <- as.numeric(ypsp$V1525) # Republicans FT, 1982
ypsp$dem97 <- as.numeric(ypsp$V5604) # Democrats FT, 1997
ypsp$rep97 <- as.numeric(ypsp$V5609) # Republicans FT, 1997
ypsp$polariz73 <- ifelse(ypsp$pid73<4,ypsp$dem73-ypsp$rep73,ifelse(ypsp$pid73>4,ypsp$rep73-ypsp$dem73,NA)) # Affective Polarization, 1973
ypsp$polariz82 <- ifelse(ypsp$pid82<4,ypsp$dem82-ypsp$rep82,ifelse(ypsp$pid82>4,ypsp$rep82-ypsp$dem82,NA)) # Affective Polarization, 1982
ypsp$polariz97 <- ifelse(ypsp$pid97<4,ypsp$dem97-ypsp$rep97,ifelse(ypsp$pid97>4,ypsp$rep97-ypsp$dem97,NA)) # Affective Polarization, 1997
## Create dataset of means and 95% CI's
library(plotrix)
figure3data <- data.frame(year=c(1973,1982,1997,1973,1982,1997),
                          Data=c("ANES","ANES","ANES","Jennings-Niemi Panel","Jennings-Niemi Panel","Jennings-Niemi Panel"),
                          means=c(mean(subset(anes.cumulative,anes.cumulative$year=="1972")$polariz.old,na.rm=T),mean(subset(anes.cumulative,anes.cumulative$year=="1982")$polariz.old,na.rm=T),mean(subset(anes.cumulative,anes.cumulative$year=="1996")$polariz.new,na.rm=T),
                                  mean(ypsp$polariz73,na.rm=T),mean(ypsp$polariz82,na.rm=T),mean(ypsp$polariz97,na.rm=T)),
                          lci=c(mean(subset(anes.cumulative,anes.cumulative$year=="1972")$polariz.old,na.rm=T) - (1.96*std.error(subset(anes.cumulative,anes.cumulative$year=="1972")$polariz.old,na.rm=T)),mean(subset(anes.cumulative,anes.cumulative$year=="1982")$polariz.old,na.rm=T) - (1.96*std.error(subset(anes.cumulative,anes.cumulative$year=="1982")$polariz.old,na.rm=T)),mean(subset(anes.cumulative,anes.cumulative$year=="1996")$polariz.new,na.rm=T) - (1.96*std.error(subset(anes.cumulative,anes.cumulative$year=="1996")$polariz.new,na.rm=T)),
                                mean(ypsp$polariz73,na.rm=T) - 1.96*std.error(ypsp$polariz73,na.rm=T),mean(ypsp$polariz82,na.rm=T) - 1.96*std.error(ypsp$polariz82,na.rm=T),mean(ypsp$polariz97,na.rm=T) - 1.96*std.error(ypsp$polariz97,na.rm=T)),
                          uci=c(mean(subset(anes.cumulative,anes.cumulative$year=="1972")$polariz.old,na.rm=T) + (1.96*std.error(subset(anes.cumulative,anes.cumulative$year=="1972")$polariz.old,na.rm=T)),mean(subset(anes.cumulative,anes.cumulative$year=="1982")$polariz.old,na.rm=T) + (1.96*std.error(subset(anes.cumulative,anes.cumulative$year=="1982")$polariz.old,na.rm=T)),mean(subset(anes.cumulative,anes.cumulative$year=="1996")$polariz.new,na.rm=T) + (1.96*std.error(subset(anes.cumulative,anes.cumulative$year=="1996")$polariz.new,na.rm=T)),
                                mean(ypsp$polariz73,na.rm=T) + 1.96*std.error(ypsp$polariz73,na.rm=T),mean(ypsp$polariz82,na.rm=T) + 1.96*std.error(ypsp$polariz82,na.rm=T),mean(ypsp$polariz97,na.rm=T) + 1.96*std.error(ypsp$polariz97,na.rm=T)))
## Create Figure 3
ggplot(figure3data,aes(x=year,y=means,color=Data,group=Data)) + geom_point() + geom_line() + theme_classic() + xlab("Year") + ylab("Affective Polarization") + geom_errorbar(aes(ymin=lci,ymax=uci),size=.5,width=0) + scale_color_manual(values=c("#000000", "#E69F00"))

# Partisan Strength
## Model
pidstrength.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='strong',age='agegroup',period='yeargroup.likes',cohort='cohort.likes.full',covariate=c(),family='binomial')
summary(pidstrength.all.raw$model) # Age and Period Effects
## Create Figure 6, Left Panel (Age Effects)
figure6leftdata <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                              effect=c(summary(pidstrength.all.raw$model)$coefficients[2:15,1]),
                              se=c(summary(pidstrength.all.raw$model)$coefficients[2:15,2]))
figure6leftdata$lci <- figure6leftdata$effect - (1.96*figure6leftdata$se)
figure6leftdata$uci <- figure6leftdata$effect + (1.96*figure6leftdata$se)
figure6leftdata$color <- ifelse(figure6leftdata$lci<0 & figure6leftdata$uci>0,"Gray50","Black")
figure6leftdata$group <- 1

figure6.a <- ggplot(figure6leftdata,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-1,1)
## Create Figure 6, Right Panel (Period Effects)
figure6rightdata <- data.frame(period=c(1955,1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                               effect=c(summary(pidstrength.all.raw$model)$coefficients[16:28,1]),
                               se=c(summary(pidstrength.all.raw$model)$coefficients[16:28,2]))
figure6rightdata$lci <- figure6rightdata$effect - (1.96*figure6rightdata$se)
figure6rightdata$uci <- figure6rightdata$effect + (1.96*figure6rightdata$se)
figure6rightdata$color <- ifelse(figure6rightdata$lci<0 & figure6rightdata$uci>0,"Gray50","Black")
figure6rightdata$group <- 1

figure6.b <- ggplot(figure6rightdata,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() +xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-1,1)
## Combine Panels
ggarrange(figure6.a,figure6.b,nrow=1)

# Ideological Sorting
## Model
ideosort.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='ideosort',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
summary(ideosort.all.raw$model) # Age and Period Effects
## Create Figure 7, Left Panel (Age Effects)
figure7leftdata <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                              effect=c(summary(ideosort.all.raw$model)$coefficients[2:15,1]),
                              se=c(summary(ideosort.all.raw$model)$coefficients[2:15,2]))
figure7leftdata$lci <- figure7leftdata$effect - (1.96*figure7leftdata$se)
figure7leftdata$uci <- figure7leftdata$effect + (1.96*figure7leftdata$se)
figure7leftdata$color <- ifelse(figure7leftdata$lci<0 & figure7leftdata$uci>0,"Gray50","Black")
figure7leftdata$group <- 1

figure7.a <- ggplot(figure7leftdata,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-0.06,0.05)
## Create Figure 7, Right Panel (Cohort Effects)
figure7rightdata <- data.frame(period=c(1980,1985,1990,1995,2000,2005,2010,2015),
                               effect=c(summary(ideosort.all.raw$model)$coefficients[16:23,1]),
                               se=c(summary(ideosort.all.raw$model)$coefficients[16:23,2]))
figure7rightdata$lci <- figure7rightdata$effect - (1.96*figure7rightdata$se)
figure7rightdata$uci <- figure7rightdata$effect + (1.96*figure7rightdata$se)
figure7rightdata$color <- ifelse(figure7rightdata$lci<0 & figure7rightdata$uci>0,"Gray50","Black")
figure7rightdata$group <- 1

figure7.b <- ggplot(figure7rightdata,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-0.06,0.05)
## Combine Panels
ggarrange(figure7.a,figure7.b,nrow=1)

# Social Sorting
## Model
socialsort.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='socialsort',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
summary(socialsort.all.raw$model) # Age and Period Effects
## Create Figure 8, Left Panel (Age Effects)
figure8leftdata <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                              effect=c(summary(socialsort.all.raw$model)$coefficients[2:15,1]),
                              se=c(summary(socialsort.all.raw$model)$coefficients[2:15,2]))
figure8leftdata$lci <- figure8leftdata$effect - (1.96*figure8leftdata$se)
figure8leftdata$uci <- figure8leftdata$effect + (1.96*figure8leftdata$se)
figure8leftdata$color <- ifelse(figure8leftdata$lci<0 & figure8leftdata$uci>0,"Gray50","Black")
figure8leftdata$group <- 1

figure8.a <- ggplot(figure8leftdata,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-0.04,0.05)
## Create Figure 8, Right Panel (Period Effects)
figure8rightdata <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                               effect=c(summary(socialsort.all.raw$model)$coefficients[16:24,1]),
                               se=c(summary(socialsort.all.raw$model)$coefficients[16:24,2]))
figure8rightdata$lci <- figure8rightdata$effect - (1.96*figure8rightdata$se)
figure8rightdata$uci <- figure8rightdata$effect + (1.96*figure8rightdata$se)
figure8rightdata$color <- ifelse(figure8rightdata$lci<0 & figure8rightdata$uci>0,"Gray50","Black")
figure8rightdata$group <- 1

figure8.b <- ggplot(figure8rightdata,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-0.04,0.05)
## Combine Panels
ggarrange(figure8.a,figure8.b,nrow=1)

# Affective Polarization among Youth vs. Rest of Electorate
## Create Data
anes.cumulative.youth <- subset(anes.cumulative,age<=7)
anes.cumulative.rest <- subset(anes.cumulative,age>=7)
polariz.youth <- data.frame(na.omit(anes.cumulative.youth %>% group_by(year) %>% summarise(means=mean(polariz.new,na.rm=T),
                                                                        lci=mean(polariz.new,na.rm=T)-(1.96*std.error(polariz.new,na.rm=T)),
                                                                        uci=mean(polariz.new,na.rm=T)+(1.96*std.error(polariz.new,na.rm=T)))))
polariz.youth$year <- c(1978,1980,1982,1984,1986,1988,1990,1992,1994,1996,1998,2000,2004,2008,2012,2016)

polariz.rest <- data.frame(na.omit(anes.cumulative.rest %>% group_by(year) %>% summarise(means=mean(polariz.new,na.rm=T),
                                                                                lci=mean(polariz.new,na.rm=T)-(1.96*std.error(polariz.new,na.rm=T)),
                                                                                uci=mean(polariz.new,na.rm=T)+(1.96*std.error(polariz.new,na.rm=T)))))
polariz.rest$year <- c(1978,1980,1982,1984,1986,1988,1990,1992,1994,1996,1998,2000,2004,2008,2012,2016)

figure12data.young <- data.frame(year=c(1978,1980,1982,1984,1986,1988,1990,1992,1994,1996,1998,2000,2004,2008,2012,2016),
                                 mean=loess(means~year,data=polariz.youth)$fitted,
                                 lci=loess(lci~year,data=polariz.youth)$fitted,
                                 uci=loess(uci~year,data=polariz.youth)$fitted)

figure12data.rest <- data.frame(year=c(1978,1980,1982,1984,1986,1988,1990,1992,1994,1996,1998,2000,2004,2008,2012,2016),
                                 mean=loess(means~year,data=polariz.rest)$fitted,
                                 lci=loess(lci~year,data=polariz.rest)$fitted,
                                 uci=loess(uci~year,data=polariz.rest)$fitted)

## Create Figure 12
plot(figure12data.young$year,figure12data.young$mean,type="l",lwd=2,ylim=c(15,45),xlab="Year",ylab="Affective Polarization")
lines(figure12data.young$year,figure12data.young$lci,lty=2)
lines(figure12data.young$year,figure12data.young$uci,lty=2)
lines(figure12data.rest$year,figure12data.rest$mean,type="l",lwd=2,ylim=c(15,45),col="blue")
lines(figure12data.rest$year,figure12data.rest$lci,lty=2,col="blue")
lines(figure12data.rest$year,figure12data.rest$uci,lty=2,col="blue")
legend("bottomright",legend=c("Rest of Electorate","Young"),fill=c("Blue","Black"))

